#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Feb 17 11:44:57 2024

@author: jcfq2
"""

import os

jwst_dir='/Users/jcfq2/data/observations/jwst'

os.chdir(jwst_dir)
import matplotlib.colors as colours
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker

import matplotlib as mpl

import scipy as sp
import scipy.ndimage

# from astropy.visualization import make_lupton_rgb
# from astropy.table import Table
# import ch4_fiddlesticks as ch4
# import pandas as pd
# import h3ppy
import glob
# import spiceypy as spice
# from astropy.io import fits
import numpy as np
from importlib import reload
# import JWSTSolarSystemPointing as jssp
# reload(jssp)
#import jwst_uranus as jwstu

# In[0]

# colors2 = plt.cm.gist_rainbow_r(np.linspace(0.25, 1, 220))
# colors1 = plt.cm.gnuplot2(np.linspace(0, 0.23, 32))


colors2 = plt.cm.jet(np.linspace(0.15, 1, int(256*(1-0.15)) ))
colors1 = plt.cm.gist_rainbow_r(np.linspace(0, 0.21, int(256*(0.15)) ))



# combine them and build a new colormap
colors = np.vstack((colors1, colors2))
mymap = colours.LinearSegmentedColormap.from_list('my_colormap', colors)


colors2 = plt.cm.gist_heat(np.linspace(0.89, 1, 56))
colors1 = plt.cm.copper(np.linspace(0, 1, 200))

# combine them and build a new colormap
colors = np.vstack((colors1, colors2))
mycoppermap = colours.LinearSegmentedColormap.from_list('my_colormap', colors)





sigma_y = 1.1
sigma_x = 1.1

null = np.nan





# In[1]




# temperature2=np.load('saturn_saves/saturn_temperature_2x.npy')
# # pre_temperature2=np.load('saturn_saves/saturn_pretemperature_2x.npy')
# temperature_error2=np.load('saturn_saves/saturn_temperature_error_2x.npy')
# density2=np.load('saturn_saves/saturn_density_2x.npy')
# # pre_density2=np.load('saturn_saves/saturn_predensity_2x.npy')
# density_error2=np.load('saturn_saves/saturn_density_error_2x.npy')
# totalE2=np.load('saturn_saves/saturn_totalE_2x.npy')
# # intensity2=np.load('saturn_saves/saturn_intensity_2x.npy')
# # ch4_fun2=np.load('saturn_saves/saturn_ch4fun_2x.npy')
# ch4_hot2=np.load('saturn_saves/saturn_ch4hot_2x.npy')
# # 

####### NB!!!! #################################################################
# since LOS wasn't working within the fits - we now fit without los and correct values here instead!!
####################################

temperature2=np.load('saturn_saves/saturn_temperature_2z.npy')
# pre_temperature2=np.load('saturn_saves/saturn_pretemperature_2.npy')
temperature_error2=np.load('saturn_saves/saturn_temperature_error_2z.npy')
density2=np.load('saturn_saves/saturn_density_2z.npy')
# pre_density2=np.load('saturn_saves/saturn_predensity_2.npy')
density_error2=np.load('saturn_saves/saturn_density_error_2z.npy')
totalE2=np.load('saturn_saves/saturn_totalE_2z.npy')
# intensity2=np.load('saturn_saves/saturn_intensity_2.npy')
# ch4_fun2=np.load('saturn_saves/saturn_ch4fun_2.npy')
ch4_hot2=np.load('saturn_saves/saturn_ch4hot_2z.npy')
# ch4_hot2=np.load('saturn_saves/saturn_ch4hot_2x.npy')
# 



####### load the LOS values ##########
naxis_maplos=np.load('saturn_spectra_los_shell_2.npy')
naxis_mapcount=np.load('saturn_spectra_count_2.npy')

nac=naxis_mapcount[:,:,600]
nal=naxis_maplos
scale=2


nal[-90*scale:-45*scale,:]=nal[-90*scale:-45*scale,:]+nal[0:45*scale,:]
nac[-90*scale:-45*scale,:]=nac[-90*scale:-45*scale,:]+nac[0:45*scale,:]
nal[45*scale:90*scale,:]=nal[45*scale:90*scale,:]+nal[-45*scale:,:]
nac[45*scale:90*scale,:]=nac[45*scale:90*scale,:]+nac[-45*scale:,:]



nal=nal[45*scale:-45*scale,:]
nac=nac[45*scale:-45*scale,:]


loscorr2=nal/nac
loscorr2=np.nan_to_num(loscorr2)


loscorr  =np.rot90(loscorr2,3)
# loscorr[:,135*scale+1:]=loscorr2[:,0:90]







temperature=np.zeros([360*scale+1,180*scale+1])
# pre_temperature=np.zeros([360*scale+1,180*scale+1])
temperature_error=np.zeros([360*scale+1,180*scale+1])
density=np.zeros([360*scale+1,180*scale+1])
# pre_density=np.zeros([360*scale+1,180*scale+1])
density_error=np.zeros([360*scale+1,180*scale+1])
totalE=np.zeros([360*scale+1,180*scale+1])
# intensity=np.zeros([360*scale+1,180*scale+1])
# ch4_fun=np.zeros([360*scale+1,180*scale+1])
ch4_hot=np.zeros([360*scale+1,180*scale+1])

temperature[:,135*scale+1:]=temperature2[:,0:90]
temperature_error[:,135*scale+1:]=temperature_error2[:,0:90]
density[:,135*scale+1:]=density2[:,0:90]
density_error[:,135*scale+1:]=density_error2[:,0:90]
totalE[:,135*scale+1:]=totalE2[:,0:90]
ch4_hot[:,135*scale+1:]=ch4_hot2[:,0:90]

#

temperature= np.fliplr(temperature.T)
density= np.fliplr(density.T)
temperature_error= np.fliplr(temperature_error.T)
density_error= np.fliplr(density_error.T)
totalE= np.fliplr(totalE.T)
ch4_hot= np.fliplr(ch4_hot.T)
# density= np.fliplr(density.T)



temperature_diff=np.zeros_like(temperature)
temperature_median = np.median(temperature,axis=1)
for xix in range(360*scale): temperature_diff[:,xix]=temperature[:,xix]-temperature_median


density_diff=np.zeros_like(density)
density_median = np.median(density,axis=1)
for xix in range(360*scale): density_diff[:,xix]=density[:,xix]-density_median
density2 = np.zeros_like(density)

for xix in range(360*scale): density2[:,xix]=density_diff[:,xix]+density_median



totalE_diff=np.zeros_like(totalE)
totalE_median = np.median(totalE,axis=1)
for xix in range(360*scale): totalE_diff[:,xix]=totalE[:,xix]-totalE_median


# %%

import cartopy.crs as ccrs
crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))

 
# Fit the non-LTE CH4 background
# crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))
crs = ccrs.NorthPolarStereo(globe=ccrs.Globe(flattening=(0.0)))
fig = plt.figure(figsize=(12,8.05),dpi=300)


bbox_props1= dict(boxstyle="circle,pad=0.15", fc="whitesmoke", ec="silver", lw=2)

fig.subplots_adjust(hspace=0.03,wspace=0.15)


# =============================================================================
# # temperature
# =============================================================================

ax = plt.subplot(331,projection=ccrs.NorthPolarStereo(central_longitude=180))

# ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())
ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())

# tp=ax.imshow(temperature, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='hot',vmin=400,vmax=480)
tp=ax.imshow(temperature, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='hot',vmin=380,vmax=480)
# tp=ax.imshow(temperature, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='hot',vmin=350,vmax=480)

ax.text(360-315,55,'a', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

cbar=fig.colorbar(tp,location='right',aspect=4)
cbar.ax.set_ylabel('Temperature [K]',labelpad=8.5)
# cbar.remove() 

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,65,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8,c='w')
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)


ax.text(360-230,64-10,'0$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-50,63-10,'180$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="right", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-329,66.5-8,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
# ax.text(360-330,69,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='mediumvioletred')
ax.text(360-144,67.7-6,'270$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
for drift in range(4): ax.plot([45+90*drift,45+90*drift],[90,40], transform=ccrs.PlateCarree(),c='hotpink',linestyle='dotted',linewidth=2)

temperature_error_zone=(temperature_error/temperature)
temperature_error_zone[temperature_error_zone == np.nan]=0
mpl.rcParams['hatch.color'] = 'w'

ax.contourf(temperature_error_zone, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),levels=[-1e19,0.05,1e19],hatches=[None,'\/',None],colors='none')
ax.contour(temperature_error_zone, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),levels=[-1e19,0.05,1e19],colors='w')


# =============================================================================
# # density
# =============================================================================


ax = plt.subplot(332,projection=ccrs.NorthPolarStereo(central_longitude=180))

ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())
# dp=ax.imshow(density, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),vmin=0,vmax=0.6e16,cmap='Greens_r')
dp=ax.imshow(density*loscorr, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),vmin=0,vmax=0.8e16,cmap='Greens_r')

ax.text(360-315,55,'b', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

cbar=fig.colorbar(dp,location='right',aspect=4)
cbar.ax.set_ylabel('Density [m$^{-2}$]',labelpad=18.5)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,65,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)


ax.text(360-230,64-10,'0$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-50,63-10,'180$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="right", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-329,66.5-8,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
# ax.text(360-330,69,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='mediumvioletred')
ax.text(360-144,67.7-6,'270$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
for drift in range(4): ax.plot([45+90*drift,45+90*drift],[90,40], transform=ccrs.PlateCarree(),c='hotpink',linestyle='dotted',linewidth=2)

density_error_zone=(density_error/density)
density_error_zone[density_error_zone == np.nan]=0
mpl.rcParams['hatch.color'] = 'k'

ax.contourf(density_error_zone, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),levels=[-1e19,0.3333,1e19],hatches=[None,'\/',None],colors='none')
ax.contour(density_error_zone, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),levels=[-1e19,0.3333,1e19],colors='k')


# =============================================================================
# # total emission
# =============================================================================

ax = plt.subplot(333,projection=ccrs.NorthPolarStereo(central_longitude=180))
# ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())
ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())

ep=ax.imshow(totalE*loscorr, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap=mymap,vmax=2e-6)

ax.text(360-315,55,'c', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

cbar=fig.colorbar(ep,location='right',aspect=4)
cbar.ax.set_ylabel('Total emission [Wm$^{-2}$]',labelpad=12)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,65,str(longit)+'$^{\circ}$W', c='w',transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)

ax.text(360-230,64-10,'0$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-50,63-10,'180$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="right", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-329,66.5-8,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
# ax.text(360-330,69,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='mediumvioletred')
ax.text(360-144,67.7-6,'270$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
for drift in range(4): ax.plot([45+90*drift,45+90*drift],[90,40], transform=ccrs.PlateCarree(),c='hotpink',linestyle='dotted',linewidth=2)







densitydiff=density_diff
densitydiff[np.abs(densitydiff)<density_error]=0

temperaturediff=temperature_diff
temperaturediff[np.abs(temperaturediff)<temperature_error]=0


# =============================================================================
# # temperature difference
# =============================================================================


ax = plt.subplot(334,projection=ccrs.NorthPolarStereo(central_longitude=180))

ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())
tp=ax.imshow(temperaturediff, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='seismic',vmin=-100,vmax=100)
# tp=ax.imshow(temperature2, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='hot',vmin=400,vmax=480)
# tp=ax.imshow(density, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),vmin=0,vmax=3e16,cmap='Greens_r')

ax.text(360-315,55,'d', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

cbar=fig.colorbar(tp,location='right',aspect=4)
cbar.ax.set_ylabel('$\Delta_{long}$ Temperature [K]',labelpad=0)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,65,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)

ax.text(360-230,64-10,'0$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-50,63-10,'180$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="right", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-329,66.5-8,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
# ax.text(360-330,69,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='mediumvioletred')
ax.text(360-144,67.7-6,'270$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
for drift in range(4): ax.plot([45+90*drift,45+90*drift],[90,40], transform=ccrs.PlateCarree(),c='hotpink',linestyle='dotted',linewidth=2)
d_temp_pos=ax.get_position()


# =============================================================================
# # density difference
# =============================================================================


ax = plt.subplot(335,projection=ccrs.NorthPolarStereo(central_longitude=180))

ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())
tp=ax.imshow(densitydiff*loscorr, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='seismic',vmin=-0.25e16,vmax=0.25e16)

ax.text(360-315,55,'e', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

cbar=fig.colorbar(tp,location='right',aspect=4)
cbar.ax.set_ylabel('$\Delta_{long}$ Density [m$^{-2}$]',labelpad=10)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,65,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)


ax.text(360-230,64-10,'0$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-50,63-10,'180$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="right", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-329,66.5-8,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
# ax.text(360-330,69,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='mediumvioletred')
ax.text(360-144,67.7-6,'270$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
for drift in range(4): ax.plot([45+90*drift,45+90*drift],[90,40], transform=ccrs.PlateCarree(),c='hotpink',linestyle='dotted',linewidth=2)
d_den_pos=ax.get_position()



# =============================================================================
# # TE difference
# =============================================================================



ax = plt.subplot(336,projection=ccrs.NorthPolarStereo(central_longitude=180))

ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())
tp=ax.imshow(totalE_diff*loscorr, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='seismic',vmin=-1e-6,vmax=1e-6)

ax.text(360-315,55,'f', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

cbar=fig.colorbar(tp,location='right',aspect=4)
cbar.ax.set_ylabel('$\Delta_{long}$ Total emission [Wm$^{-2}$]')

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,65,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)


ax.text(360-230,64-10,'0$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-50,63-10,'180$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="right", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-329,66.5-8,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
# ax.text(360-330,69,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='mediumvioletred')
ax.text(360-144,67.7-6,'270$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
for drift in range(4): ax.plot([45+90*drift,45+90*drift],[90,40], transform=ccrs.PlateCarree(),c='hotpink',linestyle='dotted',linewidth=2)




# =============================================================================
# # temperature error
# =============================================================================

ax = plt.subplot(337,projection=ccrs.NorthPolarStereo(central_longitude=180))
ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())

# tep=ax.imshow(temperature_error/temperature*100, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='hot',vmin=0,vmax=3)
tep=ax.imshow(temperature_error, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='hot',vmin=0,vmax=25)

ax.text(360-315,55,'g', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

cbar=fig.colorbar(tep,location='right',aspect=4,shrink=0.8)
# cbar.ax.set_ylabel('Temperature error %')
cbar.ax.set_ylabel('Temperature error [K]', labelpad=35)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='w', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,65,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8,c='w')



pos = cbar.ax.get_position()
cbar.ax.set_aspect('auto')

# create a second axes instance and set the limits you need
ax2 = cbar.ax.twinx()
ax2.set_ylim([0,5.2083])

# resize the colorbar (otherwise it overlays the plot)
pos.x0 +=0.02
cbar.ax.set_position(pos)
ax2.set_position(pos)
ax2.set_yticks([0,2,4],['0%','2%','4%'])



# =============================================================================
# # density error
# =============================================================================


ax = plt.subplot(338,projection=ccrs.NorthPolarStereo(central_longitude=180))
ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())

# dep=ax.imshow(density_error/density*100, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='Greens_r',vmin=0,vmax=20)
dep=ax.imshow(density_error*loscorr, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='Greens_r',vmin=0,vmax=1.5e15)

ax.text(360-315,55,'h', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

cbar=fig.colorbar(dep,location='right',aspect=6,shrink=0.8)
# cbar.ax.set_ylabel('Density error %')
cbar.ax.set_ylabel('Density error [m$^{-2}$]', labelpad=32)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,65,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)

pos = cbar.ax.get_position()
cbar.ax.set_aspect('auto')

# create a second axes instance and set the limits you need
ax2 = cbar.ax.twinx()
ax2.set_ylim([0,25])

# resize the colorbar (otherwise it overlays the plot)
pos.x0 +=0.02
cbar.ax.set_position(pos)
ax2.set_position(pos)
ax2.set_yticks([0,10,20],['0%','10%','20%'])

ax = plt.subplot(339,projection=ccrs.NorthPolarStereo(central_longitude=180))
ax.set_extent([0, 360, 60, 90], ccrs.PlateCarree())



# =============================================================================
# # methane
# =============================================================================



mp=ax.imshow(ch4_hot*1000, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),vmin=0.5,vmax=2.5,cmap=mycoppermap)
# mp=ax.imshow(ch4_hot*loscorr, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),vmin=0.000,vmax=0.001,cmap=mycoppermap)


ax.text(360-315,55,'i', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

cbar=fig.colorbar(mp,location='right',aspect=4)
cbar.ax.set_ylabel(r'CH$_4$ hot band [mW/m$^2$/sr]',labelpad=12)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,65,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)

ax.text(360-230,64-10,'0$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-50,63-10,'180$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="right", va="center",size=8,c='hotpink', weight='bold')
ax.text(360-329,66.5-8,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
# ax.text(360-330,69,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='mediumvioletred')
ax.text(360-144,67.7-6,'270$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink', weight='bold')
for drift in range(4): ax.plot([45+90*drift,45+90*drift],[90,40], transform=ccrs.PlateCarree(),c='hotpink',linestyle='dotted',linewidth=2)


fig.savefig('asd_temper_fig2.pdf', dpi=300, bbox_inches='tight', facecolor='white', pad_inches=0) 

plt.show()

# %%

# plt.plot(temperature[temperature != 0],density[temperature != 0],',')

